author=Marcio name=Listas ===== Fisica Computacional ===== Professor: [[http://complex.if.uff.br/marcio|Marcio Argollo]] Contato: Sala A1-02, tel. 2629-5773 **Referências básicas:** * Manuais de Fortran 77: * http://www.idris.fr/data/cours/lang/fortran/f90/F77.html * http://www.stanford.edu/class/me200c/tutorial_77/ * Manuais de Fortran 90 * http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_fortran90.pdf (baixar {{:apostila_fortran90.pdf|aqui}}) * http://www.cisl.ucar.edu/tcg/consweb/Fortran90/F90Tutorial/tutorial.html * Fortran 90 para programadores de Fortran 77 * http://www.nsc.liu.se/~boein/f77to90/ * Manuais de C/C++ * http://www.cenapad.unicamp.br/servicos/treinamentos/apostilas/apostila_C.pdf (baixar {{:apostila_c.pdf|aqui}}) * http://imada.sdu.dk/~svalle/courses/dm14-2005/mirror/c/ * http://www.cs.cf.ac.uk/Dave/C/CE.html * http://www.ic.unicamp.br/~cmrubira/aacesta/cpp/cpp15.html * Livros-texto * [[http://sip.clarku.edu/|An Introduction to Computer Simulation Methods: Applications to Physical Systems]], H. Gould e J. Tobochnik. Para uma revisão de comandos básicos em Unix/Linux, utilização do gnuplot e programação em C veja o tutorial apresentado na Semana da Física 2009 [[http://complex.if.uff.br/marcio/semana2009|"Ferramentas de Computação em Física"]]. ================================================================================= ==== Listas de Exercícios ==== As listas deverão ser entregues em papel e por email, compactadas no formato tar.gz e com o nome listaN_seunome.tar.gz === Lista 1) Exercicios de revisão de programação - 19/08/2010 === - Somas não são comutativas no computador! Calcule as somas S(N)=\sum_{i=1}^{N-1} \frac{1}{i^2} e SI(N)=\sum_{i=1}^{N-1} \frac{1}{(N-i)^2} com i sendo uma variável ponto flutuante de precisão simples e coloque **em um mesmo gráfico** as diferencas \left|S(N)-S_{\infty}\right|/S_{\infty} vs N e \left|SI(N)-S_{\infty}\right|/S_{\infty} vs N para N=100,200,..,10^6 onde S_{\infty}=\frac{\pi ^2}{6}. - Escrever um programa que recebe como entrada um número inteiro em representação decimal e retorna sua expansão binária. Para tal, crie uma função que verifica se o i-ésimo bit de uma palavra (variável do programa) vale zero ou um. - Considere um número de ponto flutuante de N bits com E bits reservados para o expoente e M=N-E-1 para a mantissa. Represente TODOS os números obtidos nesta representação em um grafico de barras no gnuplot. Para cada número x, armazene um par (x,1) e plote o arquivo de dados com ''with impulses''. Reveja [[http://www.psc.edu/general/software/packages/ieee/ieee.php|aqui]] a definição de números de ponto flutuante. \\ {{:floating_point.jpg?300|}} - Imagine o empilhamento de dominós de massa unitária e comprimento 2. O primeiro dominó se encontra com o centro de massa na origem, e a cada passo n deslocamos os dominós ja empilhados de \delta_n para a direita e adicionamos outro dominó na parte inferior da pilha, também com seu centro de massa na origem. Escreva um programa que diz, para um dado \delta_n=\delta~~ \forall n, quantos dominós podem ser empilhados até que o centro de massa ultrapasse a largura do primeiro dominó e a pilha tombe. (Veja uma discussão do problema [[http://web.mat.bham.ac.uk/C.J.Sangwin/Teaching/pus/dominoes.pdf|aqui]]). **A partir deste exercício todos os cálculos envolvendo reais devem ser feito em ponto flutuante com precisão dupla** === Lista 2) Derivada numérica e integração de equações diferenciais ordinárias === == O método de Euler - 31/08/2010 == Sendo a derivada de uma função f(x) definida por \frac{df}{dx}=\lim_{h \to 0}\frac{f(x+h)-f(x)}{h} e f(x+h)=f(x)+hf^{\prime}(x)+\frac{h^2}{2}f^{\prime \prime}(x)+\ldots, o método de Euler (com difereças posteriores) corresponde ao truncamento desta expansão em ordem h. Obtemos assim que a derivada no ponto x é aproximada por \frac{df}{dx}=\frac{f(x+h)-f(x)}{h}. - Mostre que o método de Euler com diferenças centradas é exato em {\cal O}(h^2), onde h é o passo de discretização da derivada. Nesse método escrevemos \frac{df}{dx}=\frac{f(x+h)-f(x-h)}{2h} - Calcule a derivada numérica de f(x)=\frac{sin(x^2)e^{x/3}}{\sqrt{(x^2+4)}} no ponto x=3 usando o método de Euler simples com passo de discretização entre h=10^{-8} e h=1 e faça um gráfico mostrando a diferença relativa entre a derivada numérica e o valor exato da derivada para ponto flutuante com precisão simples e dupla. Entre quais valores de h a aproximação é aceitável? (Solução [[http://davinci.if.ufrgs.br/wiki/index.php/Derivada_Num%C3%A9rica|aqui]]) - Mostre que o método de Euler é instável para qualquer passo h na solução do problema de crescimento exponencial \frac{dy}{dx}=\lambda x - Integre a equação diferencial \frac{dy}{dt}-y=-\frac{1}{2}e^{t/2}\sin(5t)+5e^{t/2}\cos(5t) com y(0)=1 com o método de Euler simples para diferentes valores de passo h=0.1, h = 0.05, h = 0.01, h = 0.005 e h = 0.001 e compare cada solução em t=1,2,3,4 e 5 com a solução exata y(t)=e^t+e^{t/2}sin(5t). Voce pode encontrar essa solução exata digitando a equação diferencial no site do [[http://www.wolframalpha.com|WolframAlpha]] - Faça um gráfico com a solução para h=0.05 entre t=0 e t=5 e a solução exata. - Resolva o problema acima com o método de Euler com diferenças centradas e os mesmos valores de h. Para obter o primeiro passo use o método de Euler simples com passo de discretização dez vezes menor que o usado no resto do intervalo. * Referências * http://mathworld.wolfram.com/NumericalDifferentiation.html * http://davinci.if.ufrgs.br/wiki/index.php/M%C3%A9todo_de_Euler == Integração numérica de equações diferenciais ordinárias de sistemas conservativos: métodos simpléticos - 01/09/2010== A aproximação para a derivada, como sugerida no método de Euler, introduz erros sistemáticos que acarretam na destruição de certos invariantes da dinâmica original, como a energia total no movimento do pêndulo simples. Para reduzir os efeitos da aproximação numérica à derivada, pode-se recorrer a métodos com aproximações de ordem superior (o método de Euler é de primeira ordem no passo de discretização h). == Métodos de Runge-Kutta == Ver: http://mathworld.wolfram.com/Runge-KuttaMethod.html \\ **Atenção** \\ Para sistemas de equações diferenciais, ver: http://www.nsc.liu.se/~boein/f77to90/rk.html - Resolva numericamente o problema do pêndulo físico utilizando o método de Runge-Kutta de quarta ordem com passo h=10^{-4} e ilustre o espaço de fase do sistema. Considere l=2.0 metros, onde l é o comprimento da haste de massa desprezível à qual se prende uma massa m. Use a condição inicial \theta_0=0 (posição vertical para baixo) e velocidades iniciais v_0=0.01v_c, 0.95v_c, 1.05 v_c e 1.2v_c, onde v_c é a velocidade crítica que separa os movimentos de oscilação e rotação (Ver o artigo [[http://www.sbfisica.org.br/rbef/indice.php?vol=17&num=1|"Comportamento crítico no pêndulo simples"]] do Paulo Murilo para um estudo desta transição). Considere a evolução por dois períodos do movimento. Ilustre, no mesmo gráfico, a solução do problema na aproximação de pequenos ângulos com \theta_0=0 e v_0=0.01v_c, 0.95v_c. - Faça um gráfico contendo os primeiros 30 segundos da evolução temporal dos valores analítico e numérico da energia total do pêndulo. No caso numérico, utilize o método de Euler e o método de Runge-Kutta de segunda e quarta ordem com passos h=10^{-1} e 10^{-4}. == Métodos simpléticos == Ver [[http://jcp.aip.org/resource/1/jcpsa6/v97/i3/p1990_s1?isAuthorized=no| Reversible multiple time scale molecular dynamics]] (pdf [[http://www.columbia.edu/cu/chemistry/groups/berne/papers/jcp_97_1990_1992.pdf|aqui]]) Estes métodos procuram criar esquemas de discretização da evolução Hamiltoniana de um sistema conservativo, de modo que o volume de pontos representativos no espaço de fase seja preservado durante a evolução temporal. Estes esquemas envolvem fatorização do operador Liouvilleano. - Use o método de Euler-Cromer e o método de Euler com passo h=10^{-1} e 10^{-4} para resolver numericamente o problema do pêndulo simples na aproximação de pequenas oscilações. Ulitize a condição inicial \theta_0=5^o e v_0=0. Faça um gráfico contendo a evolução temporal da posição angular por aproximadamente 10 períodos do movimento para a solução analítica e numérica com os métodos de Euler e Euler-Cromer. Faça outro gráfico ilustrando a evolução temporal da energia no mesmo alcance temporal para as 3 soluções. === Lista 3: Caos em sistemas dissipativos: o pêndulo forçado-amortecido -14/09/2010 === Nesta sequência de problemas utilize o método de Runge-Kutta de quarta ordem com passo h=10^{-4} Considere o pêndulo simples em um meio viscoso, que oferece resistência ao seu movimento com uma força proporcional à sua velocidade F_v=-\gamma v. A equação de movimento pode ser escrita como \frac{d^2\theta}{dt^2} + \frac{1}{q}\frac{d\theta}{dt} + \sin{(\theta)}=0.\\ - Faça o gráfico da posição como função do tempo para q=0.2,2.0 e 20.0 com posição angular inicial \theta_0=\pi/2 e \omega_0=0. Note que para valores grandes de q devemos recuperar o movimento sem dissipação (conservativo) e, portanto, devemos utilizar um método simplético, como Verlet ou Euler-Cromer, para integrar a equação de movimento. Considere agora uma força periódica impulsionando o pêndulo com força F_e = A\cos{(\omega_e t)}. A equação de movimento pode ser escrita como \frac{d^2\theta}{dt^2} + \frac{1}{q}\frac{d\theta}{dt} + \sin{(\theta)}=A\cos{(\omega_E t)}. - Construa o gráfico \omega(t) \times t com q=2, ~ \omega_E = 2/3 e A=1.08, 1.10, 1.13 e 1,23, respectivamente. Observe a mudança de comportamento ao variarmos a amplitude A da força externa. Simule o pêndulo por 100 ciclos da força externa e despreze os primeiros 10 ciclos (quanto vale o período \tau da força externa neste exemplo?). - Construa o espaço de fases para cada uma das situações descritas acima com a evolução do ponto representativo - Faça um gráfico "estroboscópico", chamado seção de Poincaré, imprimindo os pontos no espaço de fase correspondentes a intervalos de tempo múltiplos do período \tau da força externa t=0,\tau,2\tau,\ldots,n\tau,\ldots. - Obtenha uma série de pontos \omega_0,\omega_1,\ldots,\omega_n,\ldots onde \omega_n corresponde ao valor da velocidade angular quando \theta vale zero pela n-ésima vez durante a evolução temporal do pêndulo. Despreze os valores obtidos para tempos menores que 10 ciclos da força externa. Construa agora outra seção de Poincaré em um gráfico \omega_n \times \omega_{n+1} para os mesmos parâmetros do problema anterior. O pêndulo descrito acima deve apresentar comportamentos distintos para diferentes parâmetros. Em particular, para \omega_E=2/3 e q=2, obtemos comportamento caótico para alguns valores da amplitude da força externa A. Verifique! ^ Amplitude da força ^ Comportamento ^ | A < 1.085 | periódico | | 1.085 < A < 1.11 | caótico | | 1.11 < A < 1.14 | periódico | | 1.14 < A < 1.22 | caótico | | A ~ 1.22 | periódico | | 1.22 < A < 1.28 | caótico | | 1.28 < A < 1.475 | periódico | | 1.475 < A < 1.485 | caótico | | 1.485 < A < 1.493 | periódico | | 1.493 < A < 1.495 | caótico | | 1.495 < A < 1.497 | periódico | === Lista 4: Caos no mapa logístico - 23/09/2010 === Leia [[http://mathworld.wolfram.com/LogisticMap.html|aqui]] sobre o problema. Considere o mapa logístico x^{\prime}=f(x)=\lambda x(1-x) - Faça os gráficos (escala lin-log) da evolução temporal de x a partir de um valor inicial arbitrário para \lambda=0.01 e \lambda=0.9. - Fazer "gráficos de escada" para os mesmos valores de parâmetro e condição incial x=0.6. Como critério de parada use a condição x-x^{\prime}<10^{-3}. - Fazer "gráficos de escada" para o caso \lambda=3.1 e condição inicial x=0.3 e para a segunda iterada f(f(x)) com mesma condição inicial. - Construir o diagrama de bifurcação do mapa logístico na região 0.0<\lambda\le 4.0. - Obter os expoentes de Lyapunov \Gamma(\lambda) do mapa logístico na região 0.0<\lambda\le 4.0. O mapa logístico com \lambda=4.0 é ergódico, o que significa que podemos substituir médias temporais por médias sobre configurações de um ensemble e, nesse caso, \Gamma(\lambda)=\lim_{n\to \infty}\frac{1}{n}\sum_{i=1}^n \log{\left|f^{\prime}(x_i)\right|}=\int \log{\left|f^{\prime}(x)\right|}\rho(x) dx onde \rho(x) é a fração de vezes que a órbita se encontra entre x e x+dx. - Calcule, subdividindo o intervalo [0,1] em 100 subintervalos, o histograma \rho(x) para as primeiras 100000 iterações do mapa logístico com \lambda=4.0 e encontre \Gamma(\lambda). Esse caso particular admite a solução analítica \rho(x)=\frac{1}{\pi \sqrt{x(1-x)}}. Compare os resultados mostrando-os em um mesmo gráfico. - Obter o diagrama de bifurcação para os mapas **(a)** f(x)=r\sin{\left(\pi x\right)} **(b)** f(x)=\left(27/4\right) r x^2(1-x) ambos na região 0.75 === Lista 5: Análise de séries temporais - 19/10/2010 === - Voce encontra {{:pi_10kdigits.data.gz|neste arquivo}} a sequência dos primeiros 10 mil dígitos de \pi. * Calcular a frequência P(i) de ocorrência de cada dígito i na sequência. * Calcular a frequência P(i,j) de ocorrência consecutiva de cada par de dígitos (i,j) na primeira metade da sequência. Utilize essa informação para adivinhar, na segunda metade da sequência, qual o dígito j seguinte a um dado dígito i da sequência. Qual a sua frequência de acertos? Como ela se compara com a probabilidade medida na primeira metade? - A função de correlação c(\tau)=\frac{\left(2/T\sum_{t=0}^{T/2}x(t)x(t+\tau)\right ) - \left(2/T\sum_{t=0}^{T/2}x(t)\right )\left(2/T\sum_{t=0}^{T/2}x(t+\tau)\right )}{ \left(2/T\sum_{t=0}^{T/2}x(t)^2\right ) - \left(2/T\sum_{t=0}^{T/2}x(t)\right )^2} é uma das ferramentas tradicionais de análise de sinais. Ela mede o grau de correlação entre o valor do sinal medido em um instante e outro sinal medido em um instante de tempo subsequente. Note que, em sinais não estacionários, a função de correlação pode depender do instante inicial da medida c(t,\tau). Não trataremos desses casos aqui. - Calcule a função de correlação para uma série de 10000 pontos obtidos com a evolução do mapa logístico com \lambda=3.569945672 (ponto de acumulação, infinitas órbitas periódicas instáveis emergem neste ponto), \lambda=3.7, \lambda=1+\sqrt{2} (transição caos-ordem em uma janela de período 3) e \lambda=4.0 (ergodicidade comprovada para esse valor de parâmetro). === Lista 6: Automata celulares - 11/11/2010 === Para ler sobre automata celulares:\\ http://mathworld.wolfram.com/ElementaryCellularAutomaton.html\\ http://classes.yale.edu/fractals/CA/welcome.html\\ - Desenvolva um código para ilustrar a evolução da regra 90 utilizando somente operações lógicas e um bit para cada sítio do sistema, implementando um sistema de 32 ou 64 bits (dependendo do tamanho da palavra no computador utilizado). Implemente condições de contorno periódicas e livres e condição inicial com um bit central 1 e todos os outros 0. - Desenvolva um código para receber uma dada regra de Wolfram como parâmetro de entrada e retornar a evolução temporal de uma condição inicial por T passos de tempo em um sistema com N sítios. Utilize uma palavra de computador por sítio. Mostre graficamente as evoluções temporais para as regras 30, 60, 90, 102, 126 e 150 com N=T=200. Utilize como condição inicial um único sítio (central) com estado s=1 e todos os outros com s=0 e implemente condições de contorno periódicas (s[N]=s[0]) e livres (s[N]=s[-1]=0); - Para a regra 30, escolha um sítio e acompanhe sua evolução temporal 0,1,0,0,1,.... \\ * Calcule o número médio de bits 1 ao longo da evolução com N=10000 e T=100, 1000, 10000 e 100000.\\ * Calcule a função de correlação C(\tau)\\ * Agrupe os bits na evolução em blocos de 8 bits. Transforme o bloco em um inteiro **com sinal**. Calcule o valor médio do número obtido. Calcule a função de correlação para a série gerada com tal procedimento. Faça um gráfico com a evolução temporal dos números de 8 bits. Para lembrar da representação de inteiros com sinal leia:\\ [[http://www.cs.cornell.edu/~tomf/notes/cps104/twoscomp.html|]]\\ [[http://en.wikipedia.org/wiki/Two's_complement|]] === Lista 7: Números pseudo-aleatórios e Monte-Carlo- 23/11/2010 === Números gerados com o computador, uma máquina de Turing determinística, não pode gerar aleatoriedade com seus elementos lógicos. Podemos, no máximo, gerar sequências de números que passem por determinados [[http://www.stat.fsu.edu/pub/diehard/|testes estatísticos de aleatoriedade]]. A estes números chamamos //pseudo-aleatórios// e sua obtenção depende crucialmente do aumento de entropia gerado pela perda de bits (informação) nas operações lógicas - para os mais interessados: http://www.springerlink.com/content/jn7x3365386phn46/ Existem geradores de números verdadeiramente aleatórios, baseados na estocasticidade de processos físicos como emissão de fótons em processos radiativos ou espalhamento em divisores de feixes (veja [[http://qrbg.irb.hr/|aqui]] um exemplo). Esses geradores são usados para operações críticas, como criptografia e bingo. Uma boa fonte de informações a respeito de números pseudo-aleatórios é a biblioteca GSL, que possui, dentre diversas outras aplicações científicas como cálculo de autovalores e solução de sistemas de equações lineares, rotinas que geram números pseudo-aleatórios com diversas receitas. http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Generation.html - Usando o gerador congruencial x=(x*A)mod(B) gere pontos em um espaço 2D e visualize seu resultado. Utilize A=16807 e 65539 e B=2147483648. Estes números podem ser gerados sem a operação mod com números inteiros de 32 bits **com sinal.** Visualize seus resultados normalizando o número aleatório para gerar pontos com alcance [-1...1] - Usando o mesmo gerador com A=16807 e 65539, crie pontos em um espaço 3D. Visualize seus resultados e, girando a figura no gnuplot, tente visualizar a regularidade nos pontos gerados com 65539. Em 2D, visualize x versus 9x-6y+z onde x,y,z são a trinca de números aleatórios correspondente a um ponto no espaço 3D. Veja como o gerador 65539 é ruim para amostrar pontos no espaço. - Gere 10^4 pares de pontos (x,y) aleatoriamente distribuídos em 2D (com x e y entre -1 e 1). Calcule a fração de números gerados que satisfazem x^2+y^2 <=1 . Essa razão deve ser aproximadamente a razao entre a área de um círculo de raio unitário e um quadrado de lado 2. Esta razão vale \pi/4. Multiplique sua fração de números por 4 e voce terá calculado pi pelo Método de Monte Carlo. Faça um gráfico mostrando a convergência do valor estimado com o método de Monte Carlo para o valor aproximado de 3.14159265358979323846 (por exemplo, mostre no gráfico a diferença entre o valor aproximado e o valor estimado) utilizando o gerador com A=16807 e compare o resultado obtido com o gerador A=65539. Quantos números aleatórios devemos gerar para se obter 10\% de precisão na sua estimativa com A=16807? - Ruína do jogador: Considere dois jogadores, o primeiro com n1 moedas e o segundo com n2 moedas em um jogo de cara ou coroa. Cada vez que um jogador perde, ele entrega uma moeda para o outro jogador e o jogo se repete até um dos jogadores perder todo o dinheiro. Realize esse jogo um milhão de vezes e calcule a fração de vezes que cada jogador perdeu todo o dinheiro. Esse valor deve ser interpretado como a probabilidade P_i de cada jogador iir à falência. Compare seu resultado com os valores teóricos para as P_1=\frac{n_2}{n_1+n_2} e P_2=\frac{n_1}{n_1+n_2}. Para ler mais sobre o assunto [[http://mathworld.wolfram.com/GamblersRuin.html|clique aqui]]